#include <stdlib.h>
#include <math.h>

int sample_weights(int num_wts, double * p_wts)
{
  int i;
  double total;
  double target;
  
  total = 0;
  
  for (i=0; i<num_wts; i++)
    total += p_wts[i];

  if (0 == total)
    return rand() % num_wts;
  
  target = total * rand() / RAND_MAX;
  
  total = 0;
  
  for (i=0; i<num_wts; i++)
  {
    total += p_wts[i];
    if (total > target)
      return i;
  }

  return num_wts-1;
}

int sample_log_weights(int num_wts, double * p_log_wts)
{
  int i;
  double max_log_wt;
  
  /* compute the maximum log weight */
  max_log_wt = -HUGE_VAL;
  for (i=0; i<num_wts; i++)
    if (p_log_wts[i] > max_log_wt)
      max_log_wt = p_log_wts[i];

  if (isinf(max_log_wt))
    return rand() % num_wts;
  
  /* rescale the weights and exponentiate */
  for (i=0; i<num_wts; i++)
    p_log_wts[i] = exp(p_log_wts[i] - max_log_wt);
  
  /* now, sample from these regular weights */
  return sample_weights(num_wts, p_log_wts);
}

